Single-cell multi-ome regression models identify functional and disease-associated enhancers and enable chromatin potential analysis.

We present a gene-level regulatory model, single-cell ATAC + RNA linking (SCARlink), which predicts single-cell gene expression and links enhancers to target genes using multi-ome (scRNA-seq and scATAC-seq co-assay) sequencing data. The approach uses regularized Poisson regression on tile-level accessibility data to jointly model all regulatory effects at a gene locus, avoiding the limitations of pairwise gene-peak correlations and dependence on peak calling. SCARlink outperformed existing gene scoring methods for imputing gene expression from chromatin accessibility across high-coverage multi-ome datasets while giving comparable to improved performance on low-coverage datasets. Shapley value analysis on trained models identified cell-type-specific gene enhancers that are validated by promoter capture Hi-C and are 11× to 15× and 5× to 12× enriched in fine-mapped eQTLs and fine-mapped genome-wide association study (GWAS) variants, respectively. We further show that SCARlink-predicted and observed gene expression vectors provide a robust way to compute a chromatin potential vector field to enable developmental trajectory analysis.

Spatiotemporal kinetics of CAF-1-dependent chromatin maturation ensures transcription fidelity during S-phase.

Proper maintenance of epigenetic information after replication is dependent on the rapid assembly and maturation of chromatin. Chromatin Assembly Complex 1 (CAF-1) is a conserved histone chaperone that deposits (H3-H4)2 tetramers as part of the replication-dependent chromatin assembly process. Loss of CAF-1 leads to a delay in chromatin maturation, albeit with minimal impact on steady-state chromatin structure. However, the mechanisms by which CAF-1 mediates the deposition of (H3-H4)2 tetramers and the phenotypic consequences of CAF-1-associated assembly defects are not well understood. We used nascent chromatin occupancy profiling to track the spatiotemporal kinetics of chromatin maturation in both wild-type (WT) and CAF-1 mutant yeast cells. Our results show that loss of CAF-1 leads to a heterogeneous rate of nucleosome assembly, with some nucleosomes maturing at near WT kinetics and others showing significantly slower maturation kinetics. The slow-to-mature nucleosomes are enriched in intergenic and poorly transcribed regions, suggesting that transcription-dependent assembly mechanisms can reset the slow-to-mature nucleosomes following replication. Nucleosomes with slow maturation kinetics are also associated with poly(dA:dT) sequences, which implies that CAF-1 deposits histones in a manner that counteracts resistance from the inflexible DNA sequence, promoting the formation of histone octamers as well as ordered nucleosome arrays. In addition, we show that the delay in chromatin maturation is accompanied by a transient and S-phase-specific loss of gene silencing and transcriptional regulation, revealing that the DNA replication program can directly shape the chromatin landscape and modulate gene expression through the process of chromatin maturation.

Spatiotemporal kinetics of CAF-1-dependent chromatin maturation ensures transcription fidelity during S-phase.

Proper maintenance of epigenetic information after replication is dependent on the rapid assembly and maturation of chromatin. Chromatin Assembly Complex 1 (CAF-1) is a conserved histone chaperone that deposits (H3-H4)2 tetramers as part of the replication-dependent chromatin assembly process. Loss of CAF-1 leads to a delay in chromatin maturation, albeit with minimal impact on steady-state chromatin structure. However, the mechanisms by which CAF-1 mediates the deposition of (H3-H4)2 tetramers and the phenotypic consequences of CAF-1-associated assembly defects are not well understood. We used nascent chromatin occupancy profiling to track the spatiotemporal kinetics of chromatin maturation in both wild-type (WT) and CAF-1 mutant yeast cells. Our results show that loss of CAF-1 leads to a heterogeneous rate of nucleosome assembly, with some nucleosomes maturing at near WT kinetics and others exhibiting significantly slower maturation kinetics. The slow-to-mature nucleosomes are enriched in intergenic and poorly transcribed regions, suggesting that transcription-dependent assembly mechanisms can reset the slow-to-mature nucleosomes following replication. Nucleosomes with slow maturation kinetics are also associated with poly(dA:dT) sequences, which implies that CAF-1 deposits histones in a manner that counteracts resistance from the inflexible DNA sequence, promoting the formation of histone octamers as well as ordered nucleosome arrays. In addition, we demonstrate that the delay in chromatin maturation is accompanied by a transient and S-phase specific loss of gene silencing and transcriptional regulation, revealing that the DNA replication program can directly shape the chromatin landscape and modulate gene expression through the process of chromatin maturation.

Profiling the quantitative occupancy of myriad transcription factors across conditions by modeling chromatin accessibility data.

Over a thousand different transcription factors (TFs) bind with varying occupancy across the human genome. Chromatin immunoprecipitation (ChIP) can assay occupancy genome-wide, but only one TF at a time, limiting our ability to comprehensively observe the TF occupancy landscape, let alone quantify how it changes across conditions. We developed TF occupancy profiler (TOP), a Bayesian hierarchical regression framework, to profile genome-wide quantitative occupancy of numerous TFs using data from a single chromatin accessibility experiment (DNase- or ATAC-seq). TOP is supervised, and its hierarchical structure allows it to predict the occupancy of any sequence-specific TF, even those never assayed with ChIP. We used TOP to profile the quantitative occupancy of hundreds of sequence-specific TFs at sites throughout the genome and examined how their occupancies changed in multiple contexts: in approximately 200 human cell types, through 12 h of exposure to different hormones, and across the genetic backgrounds of 70 individuals. TOP enables cost-effective exploration of quantitative changes in the landscape of TF binding.

Cell-Cycle-Dependent Chromatin Dynamics at Replication Origins.

Origins of DNA replication are specified by the ordered recruitment of replication factors in a cell-cycle-dependent manner. The assembly of the pre-replicative complex in G1 and the pre-initiation complex prior to activation in S phase are well characterized; however, the interplay between the assembly of these complexes and the local chromatin environment is less well understood. To investigate the dynamic changes in chromatin organization at and surrounding replication origins, we used micrococcal nuclease (MNase) to generate genome-wide chromatin occupancy profiles of nucleosomes, transcription factors, and replication proteins through consecutive cell cycles in Saccharomyces cerevisiae. During each G1 phase of two consecutive cell cycles, we observed the downstream repositioning of the origin-proximal +1 nucleosome and an increase in protected DNA fragments spanning the ARS consensus sequence (ACS) indicative of pre-RC assembly. We also found that the strongest correlation between chromatin occupancy at the ACS and origin efficiency occurred in early S phase, consistent with the rate-limiting formation of the Cdc45-Mcm2-7-GINS (CMG) complex being a determinant of origin activity. Finally, we observed nucleosome disruption and disorganization emanating from replication origins and traveling with the elongating replication forks across the genome in S phase, likely reflecting the disassembly and assembly of chromatin ahead of and behind the replication fork, respectively. These results provide insights into cell-cycle-regulated chromatin dynamics and how they relate to the regulation of origin activity.

RoboCOP: jointly computing chromatin occupancy profiles for numerous factors from chromatin accessibility data.

Chromatin is a tightly packaged structure of DNA and protein within the nucleus of a cell. The arrangement of different protein complexes along the DNA modulates and is modulated by gene expression. Measuring the binding locations and occupancy levels of different transcription factors (TFs) and nucleosomes is therefore crucial to understanding gene regulation. Antibody-based methods for assaying chromatin occupancy are capable of identifying the binding sites of specific DNA binding factors, but only one factor at a time. In contrast, epigenomic accessibility data like MNase-seq, DNase-seq, and ATAC-seq provide insight into the chromatin landscape of all factors bound along the genome, but with little insight into the identities of those factors. Here, we present RoboCOP, a multivariate state space model that integrates chromatin accessibility data with nucleotide sequence to jointly compute genome-wide probabilistic scores of nucleosome and TF occupancy, for hundreds of different factors. We apply RoboCOP to MNase-seq and ATAC-seq data to elucidate the protein-binding landscape of nucleosomes and 150 TFs across the yeast genome, and show that our model makes better predictions than existing methods. We also compute a chromatin occupancy profile of the yeast genome under cadmium stress, revealing chromatin dynamics associated with transcriptional regulation.

Linking the dynamics of chromatin occupancy and transcription with predictive models.

Though the sequence of the genome within each eukaryotic cell is essentially fixed, it exists within a complex and changing chromatin state. This state is determined, in part, by the dynamic binding of proteins to the DNA. These proteins-including histones, transcription factors (TFs), and polymerases-interact with one another, the genome, and other molecules to allow the chromatin to adopt one of exceedingly many possible configurations. Understanding how changing chromatin configurations associate with transcription remains a fundamental research problem. We sought to characterize at high spatiotemporal resolution the dynamic interplay between transcription and chromatin in response to cadmium stress. Whereas gene regulatory responses to environmental stress in yeast have been studied, how the chromatin state changes and how those changes connect to gene regulation remain unexplored. By combining MNase-seq and RNA-seq data, we found chromatin signatures of transcriptional activation and repression involving both nucleosomal and TF-sized DNA-binding factors. Using these signatures, we identified associations between chromatin dynamics and transcriptional regulation, not only for known cadmium response genes, but across the entire genome, including antisense transcripts. Those associations allowed us to develop generalizable models that predict dynamic transcriptional responses on the basis of dynamic chromatin signatures.

Local nucleosome dynamics and eviction following a double-strand break are reversible by NHEJ-mediated repair in the absence of DNA replication.

We interrogated at nucleotide resolution the spatiotemporal order of chromatin changes that occur immediately following a site-specific double-strand break (DSB) upstream of the PHO5 locus and its subsequent repair by nonhomologous end joining (NHEJ). We observed the immediate eviction of a nucleosome flanking the break and the repositioning of adjacent nucleosomes away from the break. These early chromatin events were independent of the end-processing Mre11-Rad50-Xrs2 (MRX) complex and preceded the MRX-dependent broad eviction of histones and DNA end-resectioning that extends up to ∼8 kb away from the break. We also examined the temporal dynamics of NHEJ-mediated repair in a G1-arrested population. Concomitant with DSB repair by NHEJ, we observed the redeposition and precise repositioning of nucleosomes at their originally occupied positions. This re-establishment of the prelesion chromatin landscape suggests that a DNA replication-independent mechanism exists to preserve epigenome organization following DSB repair.

RoboCOP: Multivariate State Space Model Integrating Epigenomic Accessibility Data to Elucidate Genome-Wide Chromatin Occupancy.

Chromatin is the tightly packaged structure of DNA and protein within the nucleus of a cell. The arrangement of different protein complexes along the DNA modulates and is modulated by gene expression. Measuring the binding locations and level of occupancy of different transcription factors (TFs) and nucleosomes is therefore crucial to understanding gene regulation. Antibody-based methods for assaying chromatin occupancy are capable of identifying the binding sites of specific DNA binding factors, but only one factor at a time. On the other hand, epigenomic accessibility data like ATAC-seq, DNase-seq, and MNase-seq provide insight into the chromatin landscape of all factors bound along the genome, but with minimal insight into the identities of those factors. Here, we present RoboCOP, a multivariate state space model that integrates chromatin information from epigenomic accessibility data with nucleotide sequence to compute genome-wide probabilistic scores of nucleosome and TF occupancy, for hundreds of different factors at once. RoboCOP can be applied to any epigenomic dataset that provides quantitative insight into chromatin accessibility in any organism, but here we apply it to MNase-seq data to elucidate the protein-binding landscape of nucleosomes and 150 TFs across the yeast genome. Using available protein-binding datasets from the literature, we show that our model more accurately predicts the binding of these factors genome-wide.

Glucocorticoid receptor recruits to enhancers and drives activation by motif-directed binding.

Glucocorticoids are potent steroid hormones that regulate immunity and metabolism by activating the transcription factor (TF) activity of glucocorticoid receptor (GR). Previous models have proposed that DNA binding motifs and sites of chromatin accessibility predetermine GR binding and activity. However, there are vast excesses of both features relative to the number of GR binding sites. Thus, these features alone are unlikely to account for the specificity of GR binding and activity. To identify genomic and epigenetic contributions to GR binding specificity and the downstream changes resultant from GR binding, we performed hundreds of genome-wide measurements of TF binding, epigenetic state, and gene expression across a 12-h time course of glucocorticoid exposure. We found that glucocorticoid treatment induces GR to bind to nearly all pre-established enhancers within minutes. However, GR binds to only a small fraction of the set of accessible sites that lack enhancer marks. Once GR is bound to enhancers, a combination of enhancer motif composition and interactions between enhancers then determines the strength and persistence of GR binding, which consequently correlates with dramatic shifts in enhancer activation. Over the course of several hours, highly coordinated changes in TF binding and histone modification occupancy occur specifically within enhancers, and these changes correlate with changes in the expression of nearby genes. Following GR binding, changes in the binding of other TFs precede changes in chromatin accessibility, suggesting that other TFs are also sensitive to genomic features beyond that of accessibility.

MATCHER: manifold alignment reveals correspondence between single cell transcriptome and epigenome dynamics.

Single cell experimental techniques reveal transcriptomic and epigenetic heterogeneity among cells, but how these are related is unclear. We present MATCHER, an approach for integrating multiple types of single cell measurements. MATCHER uses manifold alignment to infer single cell multi-omic profiles from transcriptomic and epigenetic measurements performed on different cells of the same type. Using scM&T-seq and sc-GEM data, we confirm that MATCHER accurately predicts true single cell correlations between DNA methylation and gene expression without using known cell correspondences. MATCHER also reveals new insights into the dynamic interplay between the transcriptome and epigenome in single embryonic stem cells and induced pluripotent stem cells.

Characterization of dependencies between growth and division in budding yeast.

Cell growth and division are processes vital to the proliferation and development of life. Coordination between these two processes has been recognized for decades in a variety of organisms. In the budding yeast Saccharomyces cerevisiae, this coordination or 'size control' appears as an inverse correlation between cell size and the rate of cell-cycle progression, routinely observed in G1 prior to cell division commitment. Beyond this point, cells are presumed to complete S/G2/M at similar rates and in a size-independent manner. As such, studies of dependence between growth and division have focused on G1 Moreover, in unicellular organisms, coordination between growth and division has commonly been analysed within the cycle of a single cell without accounting for correlations in growth and division characteristics between cycles of related cells. In a comprehensive analysis of three published time-lapse microscopy datasets, we analyse both intra- and inter-cycle dependencies between growth and division, revisiting assumptions about the coordination between these two processes. Interestingly, we find evidence (i) that S/G2/M durations are systematically longer in daughters than in mothers, (ii) of dependencies between S/G2/M and size at budding that echo the classical G1 dependencies, and (iii) in contrast with recent bacterial studies, of negative dependencies between size at birth and size accumulated during the cell cycle. In addition, we develop a novel hierarchical model to uncover inter-cycle dependencies, and we find evidence for such dependencies in cells growing in sugar-poor environments. Our analysis highlights the need for experimentalists and modellers to account for new sources of cell-to-cell variation in growth and division, and our model provides a formal statistical framework for the continued study of dependencies between biological processes.

Establishment of Expression in the SHORTROOT-SCARECROW Transcriptional Cascade through Opposing Activities of Both Activators and Repressors.

Tissue-specific gene expression is often thought to arise from spatially restricted transcriptional cascades. However, it is unclear how expression is established at the top of these cascades in the absence of pre-existing specificity. We generated a transcriptional network to explore how transcription factor expression is established in the Arabidopsis thaliana root ground tissue. Regulators of the SHORTROOT-SCARECROW transcriptional cascade were validated in planta. At the top of this cascade, we identified both activators and repressors of SHORTROOT. The aggregate spatial expression of these regulators is not sufficient to predict transcriptional specificity. Instead, modeling, transcriptional reporters, and synthetic promoters support a mechanism whereby expression at the top of the SHORTROOT-SCARECROW cascade is established through opposing activities of activators and repressors.

SLICER: inferring branched, nonlinear cellular trajectories from single cell RNA-seq data.

Single cell experiments provide an unprecedented opportunity to reconstruct a sequence of changes in a biological process from individual "snapshots" of cells. However, nonlinear gene expression changes, genes unrelated to the process, and the possibility of branching trajectories make this a challenging problem. We develop SLICER (Selective Locally Linear Inference of Cellular Expression Relationships) to address these challenges. SLICER can infer highly nonlinear trajectories, select genes without prior knowledge of the process, and automatically determine the location and number of branches and loops. SLICER recovers the ordering of points along simulated trajectories more accurately than existing methods. We demonstrate the effectiveness of SLICER on previously published data from mouse lung cells and neural stem cells.

Mapping nucleosome positions using DNase-seq.

Although deoxyribonuclease I (DNase I) was used to probe the structure of the nucleosome in the 1960s and 1970s, in the current high-throughput sequencing era, DNase I has mainly been used to study genomic regions devoid of nucleosomes. Here, we reveal for the first time that DNase I can be used to precisely map the (translational) positions of in vivo nucleosomes genome-wide. Specifically, exploiting a distinctive DNase I cleavage profile within nucleosome-associated DNA--including a signature 10.3 base pair oscillation that corresponds to accessibility of the minor groove as DNA winds around the nucleosome--we develop a Bayes-factor-based method that can be used to map nucleosome positions along the genome. Compared to methods that require genetically modified histones, our DNase-based approach is easily applied in any organism, which we demonstrate by producing maps in yeast and human. Compared to micrococcal nuclease (MNase)-based methods that map nucleosomes based on cuts in linker regions, we utilize DNase I cuts both outside and within nucleosomal DNA; the oscillatory nature of the DNase I cleavage profile within nucleosomal DNA enables us to identify translational positioning details not apparent in MNase digestion of linker DNA. Because the oscillatory pattern corresponds to nucleosome rotational positioning, it also reveals the rotational context of transcription factor (TF) binding sites. We show that potential binding sites within nucleosome-associated DNA are often centered preferentially on an exposed major or minor groove. This preferential localization may modulate TF interaction with nucleosome-associated DNA as TFs search for binding sites.

Genome-wide chromatin footprinting reveals changes in replication origin architecture induced by pre-RC assembly.

Start sites of DNA replication are marked by the origin recognition complex (ORC), which coordinates Mcm2-7 helicase loading to form the prereplicative complex (pre-RC). Although pre-RC assembly is well characterized in vitro, the process is poorly understood within the local chromatin environment surrounding replication origins. To reveal how the chromatin architecture modulates origin selection and activation, we "footprinted" nucleosomes, transcription factors, and replication proteins at multiple points during the Saccharomyces cerevisiae cell cycle. Our nucleotide-resolution protein occupancy profiles resolved a precise ORC-dependent footprint at 269 origins in G2. A separate class of inefficient origins exhibited protein occupancy only in G1, suggesting that stable ORC chromatin association in G2 is a determinant of origin efficiency. G1 nucleosome remodeling concomitant with pre-RC assembly expanded the origin nucleosome-free region and enhanced activation efficiency. Finally, the local chromatin environment restricts the loading of the Mcm2-7 double hexamer either upstream of or downstream from the ARS consensus sequence (ACS).

Core and region-enriched networks of behaviorally regulated genes and the singing genome.

Songbirds represent an important model organism for elucidating molecular mechanisms that link genes with complex behaviors, in part because they have discrete vocal learning circuits that have parallels with those that mediate human speech. We found that ~10% of the genes in the avian genome were regulated by singing, and we found a striking regional diversity of both basal and singing-induced programs in the four key song nuclei of the zebra finch, a vocal learning songbird. The region-enriched patterns were a result of distinct combinations of region-enriched transcription factors (TFs), their binding motifs, and presinging acetylation of histone 3 at lysine 27 (H3K27ac) enhancer activity in the regulatory regions of the associated genes. RNA interference manipulations validated the role of the calcium-response transcription factor (CaRF) in regulating genes preferentially expressed in specific song nuclei in response to singing. Thus, differential combinatorial binding of a small group of activity-regulated TFs and predefined epigenetic enhancer activity influences the anatomical diversity of behaviorally regulated gene networks.

Convergent transcriptional specializations in the brains of humans and song-learning birds.

Song-learning birds and humans share independently evolved similarities in brain pathways for vocal learning that are essential for song and speech and are not found in most other species. Comparisons of brain transcriptomes of song-learning birds and humans relative to vocal nonlearners identified convergent gene expression specializations in specific song and speech brain regions of avian vocal learners and humans. The strongest shared profiles relate bird motor and striatal song-learning nuclei, respectively, with human laryngeal motor cortex and parts of the striatum that control speech production and learning. Most of the associated genes function in motor control and brain connectivity. Thus, convergent behavior and neural connectivity for a complex trait are associated with convergent specialized expression of multiple genes.

Learning protein-DNA interaction landscapes by integrating experimental data through computational models.

MOTIVATION: Transcriptional regulation is directly enacted by the interactions between DNA and many proteins, including transcription factors (TFs), nucleosomes and polymerases. A critical step in deciphering transcriptional regulation is to infer, and eventually predict, the precise locations of these interactions, along with their strength and frequency. While recent datasets yield great insight into these interactions, individual data sources often provide only partial information regarding one aspect of the complete interaction landscape. For example, chromatin immunoprecipitation (ChIP) reveals the binding positions of a protein, but only for one protein at a time. In contrast, nucleases like MNase and DNase can be used to reveal binding positions for many different proteins at once, but cannot easily determine the identities of those proteins. Currently, few statistical frameworks jointly model these different data sources to reveal an accurate, holistic view of the in vivo protein-DNA interaction landscape. RESULTS: Here, we develop a novel statistical framework that integrates different sources of experimental information within a thermodynamic model of competitive binding to jointly learn a holistic view of the in vivo protein-DNA interaction landscape. We show that our framework learns an interaction landscape with increased accuracy, explaining multiple sets of data in accordance with thermodynamic principles of competitive DNA binding. The resulting model of genomic occupancy provides a precise mechanistic vantage point from which to explore the role of protein-DNA interactions in transcriptional regulation. AVAILABILITY AND IMPLEMENTATION: The C source code for compete and Python source code for MCMC-based inference are available at http://www.cs.duke.edu/∼amink. CONTACT: amink@cs.duke.edu SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Stability selection for regression-based models of transcription factor-DNA binding specificity.

MOTIVATION: The DNA binding specificity of a transcription factor (TF) is typically represented using a position weight matrix model, which implicitly assumes that individual bases in a TF binding site contribute independently to the binding affinity, an assumption that does not always hold. For this reason, more complex models of binding specificity have been developed. However, these models have their own caveats: they typically have a large number of parameters, which makes them hard to learn and interpret. RESULTS: We propose novel regression-based models of TF-DNA binding specificity, trained using high resolution in vitro data from custom protein-binding microarray (PBM) experiments. Our PBMs are specifically designed to cover a large number of putative DNA binding sites for the TFs of interest (yeast TFs Cbf1 and Tye7, and human TFs c-Myc, Max and Mad2) in their native genomic context. These high-throughput quantitative data are well suited for training complex models that take into account not only independent contributions from individual bases, but also contributions from di- and trinucleotides at various positions within or near the binding sites. To ensure that our models remain interpretable, we use feature selection to identify a small number of sequence features that accurately predict TF-DNA binding specificity. To further illustrate the accuracy of our regression models, we show that even in the case of paralogous TF with highly similar position weight matrices, our new models can distinguish the specificities of individual factors. Thus, our work represents an important step toward better sequence-based models of individual TF-DNA binding specificity. AVAILABILITY: Our code is available at http://genome.duke.edu/labs/gordan/ISMB2013. The PBM data used in this article are available in the Gene Expression Omnibus under accession number GSE47026.